#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/**
 * 有理 Bézier 曲线的 de Casteljau 实现（齐次坐标法）
 *
 * 参数
 * ----
 * t : double (0 ≤ t ≤ 1)
 *     曲线参数
 * ctrl_pts : double[][]，维度为 (n+1)×2
 *     控制点列表，n 为 Bézier 曲线次数
 * weights : double[]，长度为 (n+1)
 *     对应的权重（正数）
 * n_pts : int
 *     控制点数量（n+1）
 * result : double[2]
 *     输出参数，存储 t 对应的曲线点（欧氏坐标）
 *
 * 返回值
 * ----
 * int : 0 表示成功，-1 表示参数错误
 */
int rational_de_casteljau(double t, const double *ctrl_pts, const double *weights,
                          int n_pts, double *result) {
    // 输入参数校验
    if (ctrl_pts == NULL || weights == NULL || result == NULL ||
        n_pts < 1 || t < 0.0 || t > 1.0) {
        fprintf(stderr, "Invalid input parameters\n");
        return -1;
    }

    int i, r;
    int n = n_pts - 1;  // 曲线次数
    double *homo = NULL;

    // 1. 分配齐次坐标数组内存 (n_pts × 3)：每行存储 (wx, wy, w)
    homo = (double *)malloc(n_pts * 3 * sizeof(double));
    if (homo == NULL) {
        fprintf(stderr, "Memory allocation failed\n");
        return -1;
    }

    // 初始化齐次坐标：homo[i][0] = x_i * w_i, homo[i][1] = y_i * w_i, homo[i][2] = w_i
    for (i = 0; i < n_pts; i++) {
        double w = weights[i];
        double x = ctrl_pts[i * 2];    // 控制点按行存储：(x0,y0,x1,y1,...xn,yn)
        double y = ctrl_pts[i * 2 + 1];
        homo[i * 3]     = x * w;       // wx
        homo[i * 3 + 1] = y * w;       // wy
        homo[i * 3 + 2] = w;           // w
    }

    // 2. 逐层普通线性插值（对齐次向量进行）
    for (r = 1; r <= n; r++) {
        int m = n_pts - r;  // 当前层的点数
        for (i = 0; i < m; i++) {
            // 计算当前层第 i 个点：(1-t)*homo[i] + t*homo[i+1]
            homo[i * 3]     = (1 - t) * homo[i * 3]     + t * homo[(i + 1) * 3];
            homo[i * 3 + 1] = (1 - t) * homo[i * 3 + 1] + t * homo[(i + 1) * 3 + 1];
            homo[i * 3 + 2] = (1 - t) * homo[i * 3 + 2] + t * homo[(i + 1) * 3 + 2];
        }
    }

    // 3. 投影回欧氏空间：X/W, Y/W
    double W = homo[2];  // 最终齐次点的权重分量
    if (W < 1e-12 && W > -1e-12) {  // 避免除零（权重应为正数，此处做安全检查）
        fprintf(stderr, "Zero weight encountered\n");
        free(homo);
        return -1;
    }
    result[0] = homo[0] / W;  // X/W
    result[1] = homo[1] / W;  // Y/W

    // 释放内存
    free(homo);
    return 0;
}

// 示例测试代码
int main() {
    // 测试：二次有理 Bézier 曲线（n=2，3个控制点）
    double ctrl_pts[] = {0.0, 0.0,  // P0
                         1.0, 2.0,  // P1
                         3.0, 0.0}; // P2
    double weights[] = {1.0, 2.0, 1.0};  // 对应权重
    int n_pts = 3;                       // 控制点数量
    double t = 0.5;                      // 参数 t
    double point[2];                     // 存储结果

    // 调用函数
    int ret = rational_de_casteljau(t, ctrl_pts, weights, n_pts, point);
    if (ret == 0) {
        printf("t=%.2f 对应的曲线点：(%.4f, %.4f)\n", t, point[0], point[1]);
    } else {
        printf("计算失败，错误码：%d\n", ret);
    }

    return 0;
}